home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Business Shareware
/
Business Shareware.iso
/
start
/
hobby
/
aa_51
/
ojupiter.c
< prev
next >
Wrap
C/C++ Source or Header
|
1993-02-13
|
5KB
|
240 lines
/* Orbital elements and perturbations for the planet Jupiter
* using formulas given by Meeus
*/
#include "planet.h"
#include "kep.h"
/* Calculations needed for Jupiter through Neptune:
*/
int pjup()
{
nu = T/5.0 + 0.1;
P = 237.47555 + 3034.9061*T;
Q = 265.91650 + 1222.1139*T;
S = 243.51721 + 428.4677*T;
f = 5.0*Q - 2.0*P;
V = mod360(f)*DTR;
f = 2.0*P - 6.0*Q + 3.0*S;
W = mod360(f)*DTR;
f = Q - P;
ze = mod360(f)*DTR;
f = S - Q;
psi = mod360(f)*DTR;
P = mod360(P)*DTR;
Q = mod360(Q)*DTR;
S = mod360(S)*DTR;
sinV = sin(V);
cosV = cos(V);
/* Use trigonometry to speed up sin, cos of multiple angles.
* This reduces the relative accuracy, but we only need
* absolute accuracy to about 8 places.
*/
sin2V = 2.0*sinV*cosV;
cos2V = cosV*cosV - sinV*sinV;
sinW = sin(W);
cosze = cos(ze);
sinze = sin(ze);
sin2ze = 2.0*sinze*cosze;
cos2ze = cosze*cosze - sinze*sinze;
sin3ze = sinze*cos2ze + cosze*sin2ze;
cos3ze = cosze*cos2ze - sinze*sin2ze;
sin4ze = sinze*cos3ze + cosze*sin3ze;
cos4ze = cosze*cos3ze - sinze*sin3ze;
sin5ze = sinze*cos4ze + cosze*sin4ze;
cos5ze = cosze*cos4ze - sinze*sin4ze;
sinQ = sin(Q);
cosQ = cos(Q);
sin2Q = 2.0*sinQ*cosQ;
cos2Q = cosQ*cosQ - sinQ*sinQ;
sin3Q = sinQ*cos2Q + cosQ*sin2Q;
cos3Q = cosQ*cos2Q - sinQ*sin2Q;
return(0);
}
int djup();
/* Elements and perturbations for Jupiter
*/
int ojupiter(e,J)
struct orbit *e;
double J;
{
e->epoch = J;
manoms(J);
pjup();
e->M = M5;
e->a = 5.202561;
e->ecc = (( -0.0000000017*T - 0.0000004676)*T + 0.000164180)*T + 0.04833475;
#if OFDATE
e->equinox = J;
e->i = ( 0.0000039*T - 0.0056961)*T + 1.308736;
f = (( 0.00000508*T + 0.00070405)*T + 0.5994317)*T + 273.277558;
e->w = mod360(f);
f = (( -0.00000851*T + 0.00035222)*T + 1.0105300)*T + 99.443414;
e->W = mod360(f);
f = (( -0.00000165*T + 0.0003347)*T + 3036.301986)*T + 238.049257;
e->L = mod360(f);
#else
e->equinox = J2000;
e->i = ((1.27e-7*T + 2.942e-5)*T - 0.0022374)*T + 1.305288;
f = ((8.999e-6*T - 0.00021857)*T + 0.0478404)*T + 273.829584;
e->w = mod360(f);
f = (( -1.2460e-5*T + 0.00096672)*T + 0.1659357)*T + 100.287838;
e->W = mod360(f);
f = e->M + e->w + e->W;
e->L = mod360(f);
#endif
djup(e);
return(0);
}
/* Program split arbitrarily to get it through
* Microsoft C compiler:
*/
int djup(e)
struct orbit *e;
{
double p, q;
/* perturbations in mean longitude */
p = ((-0.004692*nu - 0.010281)*nu + 0.331364)*sinV
+((0.002075*nu - 0.064436)*nu + 0.003228)*cosV
-(( -0.000489*nu + 0.000275)*nu + 0.003083)*sin2V
+0.002472*sinW
+0.013619*sinze
+0.018472*sin2ze
+0.006717*sin3ze
+0.002775*sin4ze;
f = (0.007275 - 0.001253*nu)*sinze
+0.006417*sin2ze
+0.002439*sin3ze
-(0.033839 + 0.001125*nu)*cosze
-0.003767*cos2ze;
p += f*sinQ;
f = -(0.035681 + 0.001208*nu)*sinze
-0.004261*sin2ze
+0.002178
+(-0.006333 + 0.001161*nu)*cosze
-0.006675*cos2ze
-0.002664*cos3ze;
p += f*cosQ;
f = -0.002572*sinze
-0.003567*sin2ze;
p += f*sin2Q;
f = 0.002094*cosze
+0.003342*cos2ze;
p += f*cos2Q;
/* perturbations in the perihelion */
q = (0.007192 - 0.003147*nu)*sinV
+((0.000197*nu - 0.000675)*nu - 0.020428)*cosV;
f = (0.007269 + 0.000672*nu)*sinze
-0.004344
+0.034036*cosze
+0.005614*cos2ze
+0.002964*cos3ze;
q += f*sinQ;
f = 0.037761*sinze
+0.006158*sin2ze
-0.006603*cosze;
q += f*cosQ;
f = -0.005356*sinze
+0.002722*sin2ze
+0.004483*cosze
-0.002642*cos2ze;
q += f*sin2Q;
f = 0.004403*sinze
-0.002536*sin2ze
+0.005547*cosze
-0.002689*cos2ze;
q += f*cos2Q;
e->w += q;
#if DEBUG
printf( "pperih %.5e\n", q );
#endif
f = p - (q/(e->ecc));
e->M += f;
#if DEBUG
printf( "pmanom %.5e\n", f );
#endif
e->L += p;
#if DEBUG
printf( "plong %.5e\n", p );
#endif
/* perturbations in the eccentricity */
q = ((-.0000043*nu + .0000130)*nu + .0003606)*sinV;
q += (.0001289 - .0000580*nu)*cosV;
f = -0.0006764*sinze
-0.0001110*sin2ze
-0.0000224*sin3ze
-0.0000204
+( 0.0001284 + 0.0000116*nu)*cosze
+0.0000188*cos2ze;
q += f*sinQ;
f = (0.0001460 + 0.0000130*nu)*sinze
+0.0000224*sin2ze
-0.0000817
+0.0006074*cosze
+0.0000992*cos2ze
+0.0000508*cos3ze
+0.0000230*cos4ze
+0.0000108*cos5ze;
q += f*cosQ;
f = -(0.0000956 + 0.0000073*nu)*sinze
+0.0000448*sin2ze
+0.0000137*sin3ze
+( -0.0000997 + 0.0000108*nu)*cosze
+0.0000480*cos2ze
+0.0000148*cos3ze;
q += f*sin2Q;
f = ( -0.0000956 +0.0000099*nu)*sinze
+0.0000490*sin2ze
+0.0000158*sin3ze
+0.0000179
+( 0.0001024 + 0.0000075*nu)*cosze
-0.0000437*cos2ze
-0.0000132*cos3ze;
q += f*cos2Q;
e->ecc += q;
/* perturbations in the semimajor axis */
q = -0.000263*cosV
+0.000205*cosze
+0.000693*cos2ze
+0.000312*cos3ze
+0.000147*cos4ze;
f = 0.000299*sinze
+0.000181*cos2ze;
q += f*sinQ;
f = 0.000204*sin2ze
+0.000111*sin3ze
-0.000337*cosze
-0.000111*cos2ze;
q += f*cosQ;
e->a += q;
return(0);
}
int cjupiter(e)
struct orbit *e;
{
e->plat = 0.0;
return(0);
}